Threshold dynamics of a reaction–advection–diffusion schistosomiasis epidemic model with seasonality and spatial heterogeneity

Most water-borne disease models ignore the advection of water flows in order to simplify the mathematical analysis and numerical computation. However, advection can play an important role in determining the disease transmission dynamics. In this paper, we investigate the long-term dynamics of a periodic reaction–advection–diffusion schistosomiasis model and explore the joint impact of advection, seasonality and spatial heterogeneity on the transmission of the disease. We derive the basic reproduction number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {R}}_0$$\end{document}R0 and show that the disease-free periodic solution is globally attractive when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {R}}_0<1$$\end{document}R0<1 whereas there is a positive endemic periodic solution and the system is uniformly persistent in a special case when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {R}}_0>1$$\end{document}R0>1. Moreover, we find that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {R}}_0$$\end{document}R0 is a decreasing function of the advection coefficients which offers insights into why schistosomiasis is more serious in regions with slow water flows.


Impact of schistosomiasis: worldwide and in China
Schistosomiasis, also referred to as bilharzia, is a parasitic disease which threatens public health and brings significant economic concerns, infecting approximately 200 million people worldwide.According to the World Health Organization (WHO), it has been reported in 78 countries and is endemic in 51 countries (World Health Organization 2022).
Schistosomiasis is an acute and chronic disease caused by a particular genus of blood flukes, referred to as Schistosoma.There are two main types of schistosomiasis: intestinal schistosomiasis, which is mainly caused by Schistosoma mansoni (Africa and other countries) or Schistosoma japonicum (China, Indonesia, the Philippines); and urinary schistosomiasis, caused by Schistosoma haematobium (Africa, the Middle East).Many people with schistosomiasis will not have any symptoms, or will not experience any symptoms sometimes for many months.For others, acute symptoms include fever, rash, a cough, diarrhoea, muscle and joint pain, or abdominal pain.When left untreated, eggs may travel to different parts of the body causing more serious issues.This is referred to as chronic schistosomiasis.Symptoms include anaemia, cystitis, persistent coughing, wheezing and shortness of breath, and in some cases, seizures and even paralysis.When left untreated, permanent organ damage is possible (World Health Organization 2022; National Health Service 2021).In children, who comprise approximately 123 million of all cases, the disease causes deficient growth and cognitive development (Osakunor et al. 2018).For a more complete exploration of the impacts of Schistosomiasis in humans, we refer readers to Colley et al. (2014), Rinaldo et al. (2021) While direct economic impacts of schistosomiasis is difficult to quantify, a recent study estimates that the elimination of the disease in Burkina Faso would increase average crop yields by 7-32% (Rinaldo et al. 2021).Moreover, it is noted that "villages in proximity of large-scale dams suffer an average [agricultural] yield loss of around 20%".It is therefore important to understand how proximity to water sources effects the spread of the disease.
Despite a majority of cases occurring in South Africa (World Health Organization 2022), the disease also remains endemic in some provinces of China, including Hubei, Hunan, Jiangxi, and Anhui; moreover, the monthly data of schistosomiasis cases in Hubei, Hunan and Anhui recorded by Chinese Center for Disease Control and Prevention (China CDC) display a seasonal pattern (Li et al. 2017).Of note, the Hubei province is located in the middle of the Yangtze River where the parasite is known to occupy.Global warming has also been considered as an exacerbating factor in the prevalence and spread of the disease (De Leo et al. 2020).

Mechanism of schistosomiasis transmission
Individuals become infected through direct contact with contaminated water.Transmission can also occur when people already suffering from the disease contaminate existing freshwater sources via their excrement, in which the eggs can be found.The eggs then hatch and may come in contact with the skin of other individuals using common water sources.This means that it is not possible to catch the infection directly from someone else who already has the disease.Interestingly, the transmission of schistosomiasis involves the action of an intermediate host -in this case, Oncomelania hupensis, a freshwater snail (Colley et al. 2014).Adult parasites do not only infect humans; they may also reside within mammals such as cattle, water buffalo, and many rodent species (De Leo et al. 2020).
Contamination occurs when the urine or fecal matter of infected individuals enters the water source.The eggs, contained within the excrement, then hatch as miracidia and are able to occupy the freshwater snail and develop and multiply within its intermediate host.The parasite is then able to leave the snail as cercariae and enter the fresh water source where it can survive for a number of days (Colley et al. 2014).When susceptible humans (or other mammals) wade, swim, bathe, or even perform domestic chores within the water, the cercariae are able to penetrate the skin and infect the host (World Health Organization 2022).In the preceding weeks, the parasites mature into adult worms, living in the blood vessels of the body.Here, the worms are then able to produce eggs.These eggs make their way to the bladder and/or intestine, at which point they are deposited through excrement, and the cycle may continue.Transmission can also occur if unfiltered water is taken from a contaminated body of water to be used for showers or hand washing.When left untreated, the worms are able to continue laying eggs for up to several years, or in rare cases, upwards of 40 years (Colley et al. 2014).

Prevention and control measures
Despite the debilitating impact of schistosomiasis on both individuals and communities, effective treatment and prevention measures exist and are quite well-known.The most basic means of prevention is to avoid contact with contaminated water, though this is not often an adequate solution when alternative water sources are unavailable.While no vaccine exists (Tebeje et al. 2016), other prevention measures include improving access to safe water, improved sanitation measures, educating the public on effective hygienic practices, as well as direct control of existing snail populations.For individuals already infected, Praziquantel tablets are an effective form of treatment to all species of Schitosoma (Zwang and Olliaro 2017; Centers for Disease Control and Prevention 2020), though the timing of application is important since it is most effective against adult worms (Centers for Disease Control and Prevention 2020).Some argue that new or additional drugs should be developed and utilized to best treat at risk populations (Bergquist et al. 2017).Unfortunately, an approved formulation for preschool aged children (1-4 years old) does not currently exist (World Health Organization 2022; Zwang and Olliaro 2017) and is noted as a practical limitation in the routine treatment of young children (Zwang and Olliaro 2017).
The primary method of prevention and reduction is the mass administration of drugs like Praziquantel.In high risk communities (≥ 50% prevalence in school-aged children), it is recommended to treat all school-aged children once per year along with treating adults considered at risk; for moderate risk communities (10-50% prevalence P. Wu et al. in school-aged children), treat all school-aged children once every two years along with treatment for adults in special groups only; for low risk communities (<10% prevalence), treat all school-aged children twice during their primary schooling age (Inobaya et al. 2014).This intervention is relatively low cost and yields substantial positive impact on affected communities (Inobaya et al. 2014).However, despite the low cost, the magnitude of infected individuals can make this strategy difficult to sustain.In fact, prevalence may reach equivalent levels as before the intervention within 18-24 months after ceasing treatment (Gray et al. 2010).Moreover, issues with compliance in some countries has been known to reduce the benefits of this strategy (Tallo et al. 2008).Drug tolerance has also been raised as a concern when the drug is administered repeatedly (Humphries et al. 2012).
Due to these difficulties, it is important to implement multiple measures to control schistosomiasis prevalence.One such alternative measure is to control the population of the intermediate host.A chemical compound can be used to reduce the snail population, commonly referred to as molluscicides, with varying composition and efficacy (Inobaya et al. 2014).In order for these efforts to be effective, repeated treatments are necessary, making it a time-consuming and expensive alternative, not to mention possible unintentional effects on other local species.
For other approaches, such as improved health education and behavioural changes, reasonable alternatives must be provided in order for the change to occur.For example, a water recreation area was built in Ghana to discourage youth from going into the local river, resulting in significant decrease in new cases of Schistosomiasis (Kosinski et al. 2012).Another example includes providing appropriate vessels to contain excrement from fishing boats in China in order to prevent depositing directly into the water (Wang et al. 2009).Together, these efforts result in a decrease in mortality, prevalence, and severity of infection (Inobaya et al. 2014).

Existing modelling efforts
Since the earliest schistosomiasis modelling works by Macdonald (1965) and Hairston (1965), numerous mathematical frameworks have been proposed to investigate the transmission dynamics of schistosomiasis (see, e.g., Woolhouse (1991Woolhouse ( , 1992) ) and the references therein for a systematic review of modelling efforts up to the early 1990's; see, e.g., Anderson et al. (2016) and the references therein for a more recent historical review of some existing modelling efforts).We discuss in more detail some notable efforts related to the work carried out in the current paper.
In 2002, Feng et al. (2002) formulated a schistosomiasis model with density dependence and infection age in snail population, and followed up with an estimation of parameters in 2004 (Feng et al. 2004).In 2011, Gao et al. (2011) studied a control problem for schistosomiasis transmission, considering in detail the influence of various control mechanisms on the stability of the endemic equilibrium.In 2016, to capture the effect of delays due to incubation or maturation times, Ding et al. (2016) consider a 6-dimensional system with five time delays and study the stability of the diseasefree equilibrium; Ding, Liu et al. follow up with a study of the endemic equilibrium in 2019 (Ding et al. 2019).In order to incorporate the influence of seasonality on disease transmission, Li et al. (2017) proposed a periodic model of schistosomiasis in 2017, analysing the model dynamics and performing a sensitivity analysis of the basic reproduction number with respect to some key parameters.More recently, in 2020 Zhang and Zhao (2020) developed a periodic model for schistosomiasis with time-dependent delays to study the transmission dynamics of schistosomiasis with a case study in Hubei, China.
In the models discussed above, many biological and ecological mechanisms are considered.However, the transmission mechanisms of the disease usually involve economic, psychological and social factors in addition to the intrinsic disease biology and ecology.For example, more than 82% of infected individuals lived in lake and marshland regions along the Yangtze River (Li et al. 2017).Hence, to study the impact of spatial heterogeneity and seasonality on the disease transmission, a periodic reaction-diffusion model including seasonal and spatial dynamical behaviors in the process of the disease was recently developed in Wang et al. (2018), Wang and Zhao (2008), Liang et al. (2019).As noted previously, schistosomiasis results from contact with cercariae in water, and so the impact of water flow on the survival of miracidia (cercariae) and their parasitic ability cannot be ignored.Indeed, stationary water sources contain the highest cercarial density, though the parasite may be limited by its poor ability to move horizontally (Walz et al. 2015;Upatham 1973).On the other hand, flowing water can facilitate the long distance transport of the parasite (Walz et al. 2015).In some cases, it is found that slow moving water (∼0.1 m/s) is beneficial to the spread of the parasite (Walz et al. 2015), though faster speeds (>0.13 m/s) may instead limit the ability of miracidia to infect the host snails (Upatham 1973).Therefore, advection of mircaida and cercariae as a key factor in schistosomiasis modeling should not be neglected.In fact, the relevant research works on reaction-advectiondiffusion model have begun to study the impact of advection on disease transmission (Gu et al. 2014(Gu et al. , 2015a, b;, b;Tian and Ruan 2018;Cheng and Zheng 2021).In 2018, Tian and Ruan (2018) studied the traveling wave of an advection-reaction-diffusion competition system with double free boundaries modeling invasion and competition of Aedes Albopictus and Aedes Aegypti mosquitoes.In 2021, Cheng and Zheng (2021) investigated the dynamics and spreading speed of a reaction-diffusion system with advection modeling West Nile virus.
Inspired by the above works, in this paper we consider the effect of advective forces in the periodic schistosomiasis model (Li et al. 2017), in addition to the effect of the diffusion of individuals (humans, snails, miracidia and cercariae) in heterogeneous environments.Consideration of heterogeneous environments has been applied in a wide variety of other reaction-diffusion settings, and can allow modellers to include differences in human activity or climatic factors such as local rain fall.Non-constant rates of diffusion extend further the generality and, ideally, the applicability of models to account for differences in movement behaviour, such as the influence of the presence of other organisms affecting the movement of the parasite.To the best of our knowledge, the threshold dynamics of the reaction-advection-diffusion schistosomiasis model has rarely been studied.Our goal is to gain a deeper understanding of the influence of heterogeneous environments, seasonal infection, and rates of advection on schistosomiasis transmission.
The remainder of this paper is organized as follows.We formulate our reactionadvection-diffusion schistosomiasis model with seasonality and spatial heterogeneity in Sect. 2. In Sect.3, we present the preliminaries including local and global existence of the solution and the ultimate boundedness of the system.In Sect.4, we derive the basic reproduction number R 0 of the system by applying the next generation operator.We investigate the threshold dynamics in terms of R 0 in Sect. 5.In Sect.6, we carry out numerical simulations to explore the influence of advection, seasonality and spatial heterogeneity on disease transmission dynamics, complimenting the theoretical analysis.Finally, we give a brief discussion in Sect.7.

Model formulation
To study the joint impact of seasonality and the diffusion and advection of miracidia and cercariae in rivers on the transmission of schistosomiasis, we extend the model given in Li et al. (2017) as follows: for (t, y) ∈ (0, ∞) × (0, H ) and the initial condition The growth dynamics are described by the forms introduced in Li et al. (2017).For the spatial components, we utilize Fick's law of diffusion for each component.For susceptible and infected humans and snails, and for recovered humans, we assume the flux J for each population is proportional to the gradient of each population: where d > 0 is a (possibly spatially dependent) diffusion rate corresponding to the concentration gradient of the population U .On the other hand, for the parasite populations, we assume that in addition to a diffusive flux, there is a drift downstream due to water flow proportional to the population density at some rate γ ≥ 0: For the boundary conditions, we assume that the human and snail populations satisfy a no-flux boundary condition at both ends of the river.In this case, this is equivalent to a homogeneous Neumann boundary condition at both ends of the river: On the other hand, the parasite populations are assumed to satisfy a no-flux boundary condition at the upstream end of the river, while satisfying a more general boundary condition at the downstream end of the river: (2.3)This is reasonable as the parasite populations are incapable of swimming upstream through motility alone when one considers the size of the river in relation to the size of the parasites.Then, we can describe different scenarios at the mouth of the river depending on the particular qualities of the river considered.→ ∞ capture the situation where the river stream runs into the ocean so that the downstream area environment is hostile for the survival of miracidia and cercariae.These differing scenarios will depend on the particular nature of the downstream end of the river considered; a no-flux boundary condition may only be appropriate if the downstream end of the river has a still pool resulting in a negligible loss of parasite.In reality, it is more likely that b M , b C are strictly positive.A schematic diagram of model (2.1) is shown in Fig. 1.
We now describe the components appearing in model (2.1), with a reference guide found in Table 2.1.The recruitment rate of susceptible human population is 1 and we assume that there is no immigration for infected or recovered populations.In order to incorporate the influence of seasonal variations of human activities and parasitic worm population in water as well as the spatial heterogeneity of disease transmission, we assume that the transmission rate from cercariae to humans and that from miracidia to snails are time-periodic and space-dependent functions and denote them as α 1 (t, y) and α 2 (t, y), respectively.The densities of newly infected human individuals and snails per unit time at time t and location y are given by α 1 (t, y)S 1 C C+B and α 2 (t, y)S 2 M M+K , respectively, where B and K are the corresponding half saturation constants.The natural death rate and recovery rate for humans is δ 1 and σ , respectively.Since the snail population also varies seasonally in water, we assume that both the recruitment rate ( 2 (t)) and the natural death rate (δ(t), δ 2 (t)) of snails are continuous and timeperiodic functions.Note that the death rate of infected snails δ 2 (t) may be different than that of susceptible snails δ(t), which we include as a possibility here.The population of miracidia increases when the eggs shed from infected people hatch in the water and the shedding rate is denoted by .The population of cercariae increases when they are released from snails into water and the releasing rate is denoted by ζ .In addition to the shedding from infected people and releasing from snails, we also assume logistic growth for miracidia and cercariae in water.The intrinsic growth rates of miracidia and cercariae are given by h M (t, y) and h C (t, y), respectively.The maximum capacities for miracidia and cercariae are denoted by K M and B C , respectively.The natural death rates of miracidia and cercariae are δ M and δ C , respectively.In principle, δ M and δ C could also be periodic functions of time; however, the survival time of both miracidia and cercariae are known to be significantly smaller (hours) than that of snails or humans (months or years) (Souza 2022).We therefore assume that these rates are constant in time.
The diffusion of each population is described by the term ∂ ∂ y (d j (y) ∂u ∂ y ), j = 1, 2,…, 7, respectively, where d j (y) is the location dependent diffusion rate.As noted above, our diffusive term is of Fickian form, which corresponds to neutral transition probabilities (as opposed to attractive or repulsive probabilities) (Potapov et al. 2014); other forms of diffusion may also be appropriate, but we do not explore these alternatives any further here.
Since schistosomiasis is caused by miracidia and cercariae which are swiming in rivers, the impact of water flow on their survival and movement cannot be ignored.Similar to Wang et al. (2022), we introduce the advection terms Relative rate of cercariae loss at the river's downstream end Half saturation constant of infection by miracidia (cercariae) Intrinsic growth rate of miracidia (cercariae) for miracidia and cercariae, respectively, where γ M > 0 and γ C > 0 are the advection rates.Note that positive advection rates correspond to a constant flow from upstream to downstream.

Preliminaries
From the perspective of biological application, we make the following reasonable assumptions: (A1) All model parameters are positive; and α i (t, y), i = 1, 2, are nontrivial, nonnegative, Hölder continuous functions on R × [0, H ] and ω-periodic with respect to t; (A6) 2 (t) and δ 2 (t) are nonnegative and ω-periodic with respect to t.
For convenience, we denote by dx , in which case system (2.1) can be rewritten as follows: We then have the initial and boundary conditions We first prove that the solution of system (3.1)exists locally and is positive in the positive cone of where A j (t, τ ) : X + → X + , j = 1, 2, . . ., 7, are the operators associated with and From Wang et al. (2022) (see also, e.g., Martin and Smith (1990)), we know that A j (t, τ ) is a strongly positive and compact operator on X for (t, τ ) ∈ R 2 and t > τ.Denote for ψ = (ψ 1 , ψ 2 , . . ., ψ 7 ) ∈ A + and t ≥ 0. Then system (3.1) can be rewritten as follows: where T := diag{T 1 , . . ., T 7 } and T j is given by for some α satisfying certain conditions depending on the spatial dimension and the growth rates of the system, we immediately obtain that the solution exists for all time t ∈ (0, ∞) and that there holds In our case, we will choose α = 1, which can readily be found to satisfy the necessary conditions given in the statement of Wang et al. (2022, Lemma 6.1).Our goal is to obtain for some N * (ψ) > 0.
To this end, we let = (0, H ) and set N 1 (t, y, ψ) = 3 i=1 w i (t, y, ψ).We immediately obtain the following inequality from system (3.1): Applying Gronwall's inequality we have where where Next, we set N 2 (t, y, ψ) = w 4 (t, y, ψ) + w 5 (t, y, ψ) and obtain a similar bound for N 2 .From system (3.1),we immediately obtain We then consider the following linear periodic reaction-diffusion system (3.7) From Zhang et al. (2015, Lemma 2.1), it is easy to obtain that system (3.7) has a unique ω-periodic solution S 2 (t, •) and it is globally attractive in where M 2 := M 2 (ψ) > 0; moreover we obtain lim sup t→∞ (N Next, we obtain bounds for w 6 , w 7 .Note that exp − and so one has for t ∈ (0, t ψ ).Thus, applying Gronwall's inequality once more, we obtain that there exists M 6 (ψ), M 7 (ψ) > 0 such that for t ∈ (0, t ψ ).Therefore, there exists a . Furthermore, we have shown that w j is ultimately bounded, i.e., lim sup t→∞ w i (t, y, ψ) ≤ M * for some M * > 0, uniformly for y ∈ .Finally, we define a family of maps {P(t)} t≥0 : A + → A + as P(t)ψ = w(t, •, ψ) for ψ ∈ A + .It follows from Zhao (2017) that P(t) : A + → A + is an ω-periodic semiflow.Based on the above discussions, we can verify that P(t) is ultimately bounded.Combining the compactness of P(t) with Zhao (2017, Theorem 1.1.3),we can verify that P(ω) admits a strong global attractor in A + .

The basic reproduction number of system (3.1)
In this section, we are devoted to deriving the functional expression of the basic reproduction number of system (3.1).For this purpose, we first consider the following parabolic problem It is obvious that system (4.1) has a globally attractive constant positive solution S * 1 = 1 /δ 1 ; the disease-free periodic solution is then given by E 0 = (S * 1 , 0, S 2 (t), 0, 0, 0), where S 2 (t) is the unique ω-periodic globally attractive solution of system (3.7).We first set (v 1 , v 2 , v 3 , v 4 ) = (w 2 , w 5 , w 6 , w 7 ).Then, the linearized system around E 0 is given by We denote by ) and define C ω (R, Y) as the Banach space space consisting of all ω-periodic and continuous functions equipped with the norm and Then system (4.2) can be expressed as (4.3) Obviously, for each A 5 (t, τ ), A 6 (t, τ ), A 7 (t, τ )}, t ≥ τ be the evolution operators with respect to to system dv/dt = −V (t)v.By the positivity of (t, τ ), it follows that there exists two constants > 0 and r > 0 such that || (t, τ )|| ≤ e −r (t−τ ) , t ≥ τ .For τ ≥ 0, F (τ )v(τ ) denotes the densities of newly infected individuals, infected snails, miracidia and cercariae at time τ .Then (t, τ )F (τ )v(τ ) represents the total densities of those compartments who become infected at time τ and still survive at time t.
is the density of accumulative infected individuals from time τ to time t.
If we denote the following linear operators on C ω (R, Y) by then we can verify that , which implies that L 1 and L 2 have the same spectral radius.From Zhao (2017, Chapter 3) and Bacaer and Guernaoui (2006), we can obtain that the basic reproduction number of system (3.1) is where r (L 1 ) and r (L 2 ) denote the spectral radius of L 1 and L 2 , respectively.Set Q(t) be the solution map of system (4.3) on Y and Qφ = z(t, •, φ), v = (w 2 , w 5 , w 6 , w 7 ).Then, Q(t) is the Poincaré map with respect to system (4.3).Denote by r (Q(ω)) the spectral radius of Q(ω).Based on the discussions in Wang and Zhao (2008), we have the following Lemma.
Lemma 4.1 For the basic reproduction number R 0 , one has . Then, there exists a positive ω-periodic function v * (t, y) such that v * (t, y)e λt is a solution of system (4.2).
In fact, an important research motivation of this paper is to introduce the advection term into the model to analyze the impact of water flow on disease transmission.Therefore, we give the following monotonicity proposition of basic reproduction number and the relative rates of miracidia and cercariae's loss at lower reach end of river.Proof The approach is similar to Wang et al. (2022, Proposition 3.5(ii)).We show the details for completeness.Let μ(b) be the principal eigenvalue of (4.4) with respect to b ≥ 0. We consider the monotonicity of R 0 (b) Let u (i) (t, y) be the positive eigenfunction of (4.4) with eigenvalue μ(b i ) for i = 1, 2. By the non-negativity of F , and the fact that (4.6) In fact, we can verify that u (1) (0, •) = Q (2) (ω)u (1) (0, •); otherwise, we suppose that u (1) (0, •) = lu (2) (0, •) for some l > 0, with r (Q (2) (ω)) = 1 the principle eigenvalue with respect to u (2) (0, •).

The global dynamics of system (3.1)
To study the global dynamics of system (3.1),we first give the following lemma which will be used later.

The global asymptotic stability of the disease-free periodic solution E 0
We first state the following Theorem.
Proof To complete the proof, we first introduce the following notations to be used later.

The uniform persistence of the disease
Theorem 5.3 If R 0 > 1, then there exists a constant > 0 such that for any ψ ∈ A + and at least one ψ i ≡ 0 for i = 2, 5, 6, 7, we have Furthermore, system (3.1)admits at least one positive ω-periodic solution.
Proof Consider the following system with > 0 for (t, y) ∈ [0, ∞) × (0, H ) with boundary condition For any φ ∈ Y, let be the unique solution of system (5.6) with v (t, •, φ) = φ.Let Q (ω) : Y → Y be the Poincaré map of system (5.6) and denote the spectral radius of If it is not true, then there exists (5.7) From Lemma 5.1 and (5.7), we have where v (n 2 ω, y) = e λ t v * (t, y), and v * (t, y) is a positive ω-periodic function associated with the solution v (n 2 ω, y) of system (5.6) with λ = ln r (Q (ω))/ω.Hence, we have which leads to a contradiction with the ultimately boundedness of system (3.1).

Numerical simulation
In this section we explore some numerical simulations to compliment the preceding analysis.In particular, we are interested in the influence of heterogeneity of the environment and advection speeds on the transmission of schistosomiasis.For simplicity, we choose and the initial values are chosen to be where 173 × 20 7 , 124, 30, 7.547 × 10 5 , 2.632 × 10 5 , 259162.05 × 10 7 ).
These parameter ranges for the initial data are based upon the reported data found in Li et al. (2017).For estimating diffusion rates, acquiring precise values or even a suitable order of magnitude for each population poses a significant challenge.Consequently, the chosen rates, though somewhat arbitrary, are primarily aimed at providing a more insightful representation of the analytical findings.

Setup of the finite difference method
In this subsection, we apply the finite difference method to achieve discretization of model (3.1) in one-dimensional space.This numerical method has been used in many PDE infectious disease model studies (Tian and Ruan 2018;Wang et al. 2022), and convergence and accuracy can be guaranteed.According to Table 2, we can rewrite w j (t, y) (i = 1, 2, . . ., 5) equations model (3.1) as follows u y (t, 0) = u y (t, H ) = 0, t > 0, (6.1)where u = (w 1 , w 2 , w 3 , w 4 , w 5 We divide the interval [0, H ] into n equally sized sub-intervals [y i , y i+1 ], i = 1, . . ., n + 1 with y i = y i−1 + y, y 1 = 0, y n+1 = L and y = H /n. We also denote the step size by t = T max /m, j = 1, . . ., m. Set u i j = u(t j , y i ), t + j = t j−1 + t, y i = y i−1 + y, 1 ≤ j ≤ m, 1 ≤ i ≤ n + 1.We use backward and central differences to discretize u t , u y and u yy as follows: Thus, for w j (t, y) in model (6.1), we have For the boundary points i = 1 and i = n + 1, since model (3.1) satisfies the Neumann boundary condition, it follows that u y (0) = 0, u y (t j , 0) = u 1 j −u 0 j y = 0, so we have u 0 j = u 1 j .Then, Similarly, we can rewrite w 6 (t, y) and w 7 (t, y) in model (3.1) as follows Thus, we can obtain As argued previously, since model (6.2) satisfies Neumann and no-flux boundary conditions, it follows that dv y (t, 0) = 0, and v 0 j = v 1 j .Thus, we have
Note that Lemma 4.3 provides us a method to calculate the value of R 0 numerically.To this end, set ϕ ∈ Int(Y + ) and denote by for m ≥ 1 and any μ ∈ (0, ∞).Here {W (t, s, μ) : t ≥ s} is the family of evolution operators F and V associated with the linear equation ), we know that r (W (ω, 0, μ)) = lim m→∞ g m whenever the limit exists.Then, we use that R 0 = 1/μ.Utilizing the MATLAB software, we find that R 0 = 0.7912 < 1, which implies that schistosomiasis transmission cannot be sustained (see Fig. 2), which aligns with the theoretical result of Theorem 5.2.When we choose α 1 = 8 × 10 −14 , α 2 = 1.974 × 10 −8 and γ M = 0.1, γ C = 0.2 and all other parameter values as in Fig. 2. We then find that R 0 = 2.518 > 1 and so schistosomiasis and miracidia(cercariae) prevalence oscillates in both river and human hosts.That means that schistosomiasis will persist in both river and human hosts (see Fig. 3), which aligns with the theoretical result of Theorem 5.3.

The impact of advection, spatial heterogeneity and seasonality on the schistosomiasis transmission
To study the impact of advection on the schistosomiasis spread, similar to the numerical method found in Aatif et al. (2023), we examine the influence of the advection rates γ M and γ C on the value of R 0 in Fig. 4a and b, respectively.It can be seen that R 0 decreases sharply with γ M ranging from 0 to 1 and with γ C ranging from 0 to 2, respectively.The value of R 0 goes down slowly as γ M and γ C continue to increase to 4. This indicates that a faster flow rate of the river could lower the risk of a schistosomiasis outbreak.The corresponding impact of advection rates of miracidia and cercariae on the dynamical behavior of system (3.1) is shown in Fig. 6.
We are also interested in how spatial heterogeneity of the transmission rates from cercariae to humans and from miracidia to snails affects the value of R 0 .To this end, we assume α i (t, y) = α i F i (y) • T i (t) (i = 1, 2), where T i (t) = 1 + 0.5 sin(0.25 + π t 6 ), and F i (y) = 0.85 + Fig. 3 The dynamical behavior of system (2.1) when R 0 > 1 0.8 cos(2π y).To study the influence of F(y) on R 0 , we set F i (y; k; s) = 0.85+k0.8cos(sπ y) with k ∈ [0, 1], s = 1, 2. Figure 4c and d show the changes in the spatial dependence of α i (t, y) decreases and increases the risk of the disease outbreak for s = 1 and s = 2, respectively.We set T i (t) = m(1 + 0.5 sin(0.25 + 2π t ω )) to study the impact of seasonality on the value of R 0 in Fig. 4. In Fig. 4e, the influence of seasonality is shown by comparing the value of R 0 under annual infection ω = 12, semi-annual infection ω = 6 and quarter-infection ω = 3, respectively.It is obvious that the disease risk R 0 is an increasing function with respect to m and the risk of the disease outbreak under annual infections is apparently higher than in quarter-infection and semi-infection cases.This indicates that seasonality plays a key role in the transmission of schistosomiasis.We also examine the impact of half saturation rate B of cercariae on the disease risk R 0 in Fig. 4d.As expected from the form appearing in the model, we observe that R 0 decreases as B increases (i.e., as B increases the infection rate decreases), which implies the saturation of infection function also has a certain impact on the spread of disease: this value could increase or decrease depending on the level of resistance to infection in the human population or the parasite's ability to infect the human population.

Sensitivity analysis
The main objective of this subsection is to discuss the sensitivity of R 0 .For this purpose, we apply the normalized forward sensitivity index and the central difference approximation to evaluate the derivatives (Chitnis et al. 2008;Kong et al. 2018;Wang et al. 2019) as follows: (6.3) Fig. 4 The impact of advection, spatial heterogeneity and seasonality on R 0 .The values of the other parameters are the same as in Fig. 3 Setting h = 1% of the parameter value, Equation (6.3) becomes S.I. = R 0 (1.01P) − R 0 (0.99P) 0.02R 0 (P) .
Using this formula we obtain that the sensitivity indices of R 0 to parameters γ M , γ C , λ 1 , 2 , σ , , ζ , α 1 , α 2 , δ 1 , δ 2 , δ M and δ C are, respectively, −0.6432, −0.6066, 1.033, 0.844, −0.3008, 0.2345, 0.4091, 0.6908, 0.5778, −0.3806, −0.2088, −0.233, and −0.344 (see Fig. 5).In Fig. 5, we fixed the parameter values when we calculate the sensitivity indices.It Fig. 5 The impact of advection on the dynamical behavior of the solution of system (3.1).The values of the other parameters are the same as in Fig. 3 can be seen that the parameters 1 , 2 , , ζ, α 1 , α 2 are positively correlated with R 0 , while γ M , γ C , σ, δ 1 , δ 2 , δ M , δ C are negatively correlated with R 0 .This indicates that increasing the recruitment rate of sea snails, the shedding rate of parasitic eggs and the releasing rate of cercariae from snails will increase the value of the basic reproduction number R 0 , while increasing the mortality rate of sea snails, as well as the mortality rate of larvae and cercariae, will reduce the value of the basic reproduction number.In practical terms, the removal of snails, larvae, and cercariae in rivers will effectively reduce the outbreak of diseases, which is actually one of the most effective methods commonly adopted for schistosomiasis prevention and control.Of particular interest is the relatively strong negative correlation with the advection rates γ M and γ C : given the option, it may be advisable to engage in bathing/cleaning activities in areas of the river with a moderate flow rate as opposed to areas of still water.This is a simple and cost effective approach when compared to other, more invasive avenues of reduction.

Discussion
Due to the complexity of the disease and the vast variety of environmental conditions in which it is found and contracted, the mathematical modeling and dynamical analysis of schistosomiasis plays a crucial role in furthering our understanding of its spread and evaluating potential mitigation measures.In this paper, we formulated a novel, non-autonomous reaction-advectiondiffusion model that incorporates hosts (susceptible infected and recovery individuals), snails (healthy and infected snails), miracidia and cercariae, to investigate the impact of spatial heterogeneity and seasonality on the transmission of schistosomiasis.Moreover, as a difference, we add drift terms to describe the effect of advection in water on miracidia and cercariae due to river flow, which has been ignored in previous modelling efforts.Applying the operator semigroup theory and the conception of the next generation operator, we obtain the key threshold found in many epidemic models: the basic reproduction number R 0 .We then study the asymptotic behavior of R 0 for the advection and diffusion rates.Detailed analyses of the global stability  2 of the disease-free periodic solutions and the uniform persistence of the disease are provided by using theory of dissipative and monotonic systems.
To compliment these analytical findings, we also investigate changes in R 0 and solution behaviours for some specific cases using numerical simulation.From these simulations, we have observed the following findings and their related biological meaning: (i) We conduct the numerical simulations of the dynamical behavior of system (3.1) in Figs. 2 and 3, which is in lines with the theoretical result mentioned in the corresponding Theorem.(ii) Numerical simulations indicate that seasonality and spatial heterogeneity are two key factors that impact the transmission of schistosomiasis; in Fig. 4, the sensitivity analysis results of R 0 for some parameters illustrate that spatial heterogeneity will amplify or dilute schistosomiasis transmission, and it is highly associated with parameters.Moreover, it indicates that seasonality together with the spatial heterogeneity may increase the complexity of schistosomiasis transmission.(iii) We investigate the role of the advection rates γ M and γ C of miracidia and cercariae on the basic reproduction number R 0 and dynamical behavior of schistosomiasis outbreak in Figs.4a, b and 6, respectively.We find that the advection rates γ M and γ C can reduce the value of R 0 to some extent and we also notice that the cumulative number of infected individuals decreases with convection rates γ M and γ C , which implies that schistosomiasis is generally difficult to break out in fast flowing waters.In other words, the key areas for schistosomiasis prevention and control should be the waters with gentle water flow; alternatively, a cost effective mitigation measure may be to simply encourage bathing/cleaning behaviours in regions with moderate water flow.We were not able to take into account the specific measures such as drug treatment, snail control, cercariae control, improved sanitation and health education, but we will leave such problems for our future study.
This work expands upon previous modelling efforts through inclusion of spatial effect and other environmental factors not previously considered.More precisely, the founding works (Macdonald 1965;Hairston 1965) and subsequent efforts focus primarily on a spatially homogeneous (i.e. through ordinary differential equations) setting.In cases where a partial differential equation setting is used, as in the works (Chan et al. 1995;Zhang et al. 2007), for example, it is most commonly motivated by an age-structure modelling framework.In this sense, the addi-tional "space" dimension included is not physical space as treated in the present work.Motivated by the formulation found in Li et al. (2017), we extend their model to include physical space in addition to the influence of seasonality.In our case, space is taken to be a one-dimensional, bounded domain (0, H ), with boundary conditions chosen such that the model corresponds most closely with schistosomiasis spread in a river system.Naturally, diffusion is a key factor, which includes cases of spatially dependent diffusion processes.Through its construction, we are also able to include advective forces caused by flowing water in the river.Inclusion of such environmental effects are important to further our understanding of disease spread in areas such as villages along the Yangze River region.One interesting finding is the effect of slow or fast flowing water on disease spread; our findings are consistent with case studies, where schistosomiasis is less prevalent in areas with faster flowing water when compared to areas with slower water flow.Such findings can provide further insights into less intrusive mitigation measures.For example, chemicals used to eliminate the snail population from sources of fresh water may have a negative impact on other local species.Furthermore, without continual treatment of the water source, it is likely that the snail population will inevitably return (Centers for Disease Control and Prevention 2020).
Despite the novel extensions presented here, there are numerous areas ripe for further study.As performed in the works of Chan et al. (1995), Zhang et al. (2007), it would be interesting to consider an age-structure in addition to spatial effects considered here.Our work also does not include many aspects that we know to be relevant to disease transmission: social, psychological, and economic factors should not be entirely ignored.It would be a challenging but potentially fruitful task to incorporate these additional effects to determine optimal control strategies as they relate to explicit spatial effects.More broadly, the modelling approach presented here could be applied to a wider range of scenarios which include the transmission of disease via a river system.
The parameters b M and b C respectively denote the relative loss rates of miracidia and cercariae due to water flux at the river's downstream end.In particular, b M = 0 and b C = 0 represent the no-flux boundary condition; b M = 1 and b C = 1 yield the homogeneous Neumann boundary condition; and b M → ∞ and b C → ∞ lead to the homogeneous Dirichlet boundary condition (i.e., M(t, H ) = C(t, H ) = 0).The case of b M → ∞ and b C

Fig. 1
Fig. 1 The schematic diagram of model (2.1) By the comparison principle and Teng et al. (2008, Lemma 1), there holds

Fig. 6
Fig. 6 Sensitivity of R 0 to parameters.The value of parameters used in computation of sensitivity indices are the same as in Table 2

Table 1
The biological interpretations of the parameters of system (2.1) MRelative rate of miracidia loss at the river's downstream end b C Proof We first set t ψ be the maximal existence time of the solution w.We follow the proof ofWang et al. (2022, Theorem 2.1).More precisely, byWang et al. (2022,  Lemma 6.1) if we can show that there exists a positive function M α (ψ) such that 0 = ψ.Furthermore, the ω-periodic semiflow P(t) on A + generated by system (3.1) is defined as P(t)ψ = w(t, •, ψ), and P(ω) admits a strong global attractor in A + , t ≥ 0.

Table 2
The value of the parameters of system (2.1)